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Abstract 


DRDC  -  Atlantic  Research  Centre,  an  agency  of  the  Canadian  Department  of  National 
Defence,  has  developed  a  prototype  multi-influence  underwater  surveillance  sensor 
called  the  Starfish.  The  Starfish  has  acoustic,  electric,  and  magnetic  held  sensors,  and 
onboard  signal  processing  for  target  detection  and  tracking.  The  acoustic  sensors  are 
in  the  form  of  a  seven-element  cross-dipole  hydrophone  array,  with  two  hydrophones 
along  each  orthogonal  axis  direction  and  one  at  the  array  centre. 

Acoustic  target  detection  is  performed  using  both  time-domain  and  frequency-domain 
Sentinel  algorithms.  The  Sentinel  algorithms  use  parallel  exponential  filters,  with  the 
frequency-domain  algorithm  using  parallel  filters  in  each  frequency  bin. 

A  Phase  Gradient  bearing  estimation  algorithm  is  used  to  perform  the  acoustic  target 
tracking.  This  is  a  frequency-domain  algorithm  that  estimates  the  gradient  of  the 
spatial  phase  of  the  incoming  acoustic  signal,  and  directly  computes  the  wavevector  of 
the  acoustic  arrival  in  each  frequency  bin.  A  Histogram  Tracking  procedure  is  used  to 
display  and  determine  the  track  angles  for  a  broad-band  acoustic  target.  This  paper 
will  discuss  the  details  of  the  frequency-domain  Acoustic  Sentinel  algorithm  and  the 
Phase  Gradient  bearing  estimation  algorithm  with  Histogram  Tracking. 

Significance  for  defence  and  security 


For  the  Force  ASW  project,  the  DRDC  Starfish  sensor  units  provide  underwater 
sensor  capability  for  acoustic  target  detection  and  tracking.  The  DRDC  Starfish 
units  use  a  seven-element  three-dimensional  acoustic  hydrophone  array  to  sample 
the  local  underwater  acoustic  environment,  and  have  onboard  signal  processing  to 
determine  the  presence  and  track  of  a  broad-band  acoustic  target.  The  Starfish 
autonomous  sensor  systems  with  their  detection  and  tracking  algorithms  give  the 
Canadian  Forces  surveillance  capability  to  detect  surface  and  subsurface  acoustic 
targets  in  littoral  waters  and  choke  points. 
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Resume 


Le  Centre  de  recherches  de  l’Atlantique  de  RDDC,  une  agence  du  ministere  de  la  De¬ 
fense  nationale,  a  mis  au  point  un  prototype  de  captenr  de  surveillance  sous-marine 
a  influences  multiples  appele  Starfish.  Le  Starfish  est  muni  de  capteurs  acoustiques, 
electriques  et  de  profondeur  de  champ  magnetique  en  plus  d’une  unite  embarquee 
de  traitement  du  signal  pour  la  detection  et  la  poursuite  de  ciblcs.  Les  capteurs 
acoustiques  sont  const  it  ues  d’un  reseau  hydrophonique  de  dipoles  croises  a  sept  com- 
posantes,  soit  deux  hydrophones  dans  chaque  axe  orthogonal  et  un  autre  au  centre 
du  reseau. 

Les  ciblcs  acoustiques  sont  detectees  au  rnoyen  des  algorithmes  Sentinel  dans  le  do- 
maine  temporel  et  le  domaine  frequentiel.  Ces  algorithmes  font  appel  a  des  Litres 
exponentiels  paralleles ;  dans  le  domaine  frequentiel,  l’algorithme  utilise  des  Litres 
parallcles  dans  chaque  cellule  de  frequences. 

Lin  algorithmic  d’estimation  d’azimut  du  gradient  de  phase  sert  a  la  poursuite  des 
cibles  acoustiques.  11  s’agit  d’un  algorithmic  dans  le  domaine  frequentiel  qui  permet 
d’evaluer  le  gradient  de  la  phase  spatiale  du  signal  acoustique  entrant  et  de  calculer 
directement  le  vecteur  d’onde  d’arrivee  dans  chaque  cellule  de  frequences.  Line  proce¬ 
dure  de  poursuite  par  histogramme  sert  a  determiner  et  a  afhcher  les  angles  de  route 
des  cibles  acoustiques  a  large  bande.  Le  present  article  expose  en  detail  1’algorithme 
acoustique  Sentinel  dans  le  domaine  frequentiel  et  l’algorithme  devaluation  d’azimut 
du  gradient  de  phase  avec  poursuite  par  histogramme. 


Importance  pour  la  defense  et  la  securite 


Pour  le  projet  de  guerre  anti-sous-marine  (GASM)  de  la  force,  le  Starbsh  de  RDDC 
fournit  une  capacite  sous-marine  de  detection  et  de  poursuite  de  cibles  acoustiques. 
Les  capteurs  Starbsh  de  RDDC  utilisent  un  reseau  hydrophonique  tridimensionnel  a 
sept  composasntes  pour  echantillonner  l’environnement  acoustique  sous-marin,  et  ils 
sont  dotes  d’une  unite  embarquee  de  traitement  du  signal  pour  determiner  la  presence 
d’une  cible  acoustique  a  large  bande  et  la  poursuivre.  Les  systemes  autonomes  de 
detection  Starbsh  et  leurs  algorithmes  de  detection  et  de  poursuite  procurent  aux 
Forces  arrnees  canadiennes  une  capacite  de  surveillance  permettant  de  detecter  les 
cibles  acoustiques  a  la  surface  et  sous  la  surface  des  eaux  littorales  et  des  points  de 
passage  oblige. 
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1  Introduction 


The  DRDC  -  Atlantic  Research  Centre  has  had  a  multi-year  research  and  development 
program  to  create  a  prototype  underwater  multi-influence  sensor  package  for  littoral 
surveillance  applications.  The  prototype  sensors,  shown  in  Fig.  1,  are  in  their  second 
generation  and  are  called  Starfish  sensor  nodes.  The  Starfish  nodes  are  designed  to  be 
deployed  in  up  to  500  m  water  depth.  Each  Starfish  sensor  node  contains  a  low-noise, 
high  sensitivity  tri-axial  magnetometer,  three  Underwater  Electric  Potential  sensors, 
a  seven-element  crossed-dipolc  hydrophone  array,  a  wide-band  hydrophone,  and  a 
combined  3D  gravitational  and  magnetic  orientation  sensor.  Each  Starfish  node  con¬ 
tains  a  powerful  onboard  ARM  processor  that  runs  an  embedded  version  of  the  Linux 
operating  system.  The  processor  is  powerful  enough  to  perform  real-time  acquisition, 
storage,  and  signal  processing  of  the  sensor  data.  Target  detection  and  tracking  re¬ 
sults  from  the  onboard  processing  are  relayed  through  an  attached  cable  or  a  network 
of  underwater  acoustic  modems  to  a  surface  gateway  buoy,  which  has  a  wireless  radio 
frequency  link  to  shore.  Each  Starfish  node  contains  internal  solid-state  memory  stor¬ 
age  where  sensor  data  may  also  be  logged  for  later  analysis.  The  Starfish  nodes  may 
be  deployed  with  or  without  cables  to  shore.  When  deployed  without  cables,  each 
Starfish  node  can  run  autonomously  for  more  than  five  days  of  continuous  use,  or  for 
several  weeks  in  a  wake-on-detection  mode.  This  paper  describes  the  acoustic  pro¬ 
cessing  running  onboard  the  Starfish  sensor  node.  The  acoustic  processing  proceeds 
as  follows:  first  the  acoustic  Sentinel  algorithm  determines  if  there  is  an  acoustic  tar¬ 
get  nearby,  second  the  acoustic  tracking  algorithm  using  the  Phase  Gradient  bearing 
estimation  algorithm  with  Histogram  Tracking  determines  the  target’s  bearing  and 
bearing  change  rate,  and  finally  a  compact  NILUS  protocol  message  [1]  is  generated 
and  transmitted  to  the  user  over  the  communications  network. 


2  Frequency-domain  acoustic  sentinel 


The  Starfish  Cubes  have  two  algorithms  that  determine  if  an  acoustic  target  is  present, 
the  time-domain  and  the  frequency-domain  Acoustic  Sentinel.  The  time-domain 
Acoustic  Sentinel  calculates  the  average  power  from  the  array  centre  hydrophone, 
and  passes  this  value  into  two  exponential  filters,  one  with  a  long  time  constant 
and  one  with  a  shorter  time  constant.  The  ratio  of  the  outputs  of  the  fast  to  slow 
filters  is  used  to  determine  if  a  target  is  present.  Similarly,  the  frequency-domain 
Acoustic  Sentinel  algorithm  works  on  Fast  Fourier  Transform  (FFT)  processed  data 
collected  from  the  centre  hydrophone  of  the  crossed-dipole  hydrophone  array.  Each 
channel  of  the  hydrophone  array  is  sampled  at  a  rate  of  2500  samples  per  second, 
and  is  low-pass  filtered  with  a  700  Hz  cutoff  frequency.  The  typical  frequency  band 
for  frequency-domain  Acoustic  Sentinel  processing  is  50  Hz  to  500  Hz.  The  frequency- 
domain  Sentinel  algorithm  uses  a  bank  of  exponential  digital  filters,  two  filters  for 
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each  FFT  frequency  bin  in  the  processing  bandwidth.  The  input  to  each  pair  of  fil¬ 
ters  is  the  average  power  spectrum  value  obtained  from  the  FFT  processing  for  that 
frequency  bin.  The  general  input-output  relationship  of  the  exponential  filter  is  given 
by  Eqn.  (1), 


y(n)  =  (1  —  a)x{n)  +  ay(n  —  1),  n>  0  (1) 

n 

=  (1  —  a)  ^2  oikx(n  —  k ) 
k= o 

and  where  x{n)  =  y(n)  =  0  for  n  <  0.  In  Eqn.  (1),  x{n)  and  y(n)  are  the  hlter 
input  and  output  respectively  at  sample  n.  The  parameter  a  determines  the  low-pass 
cut-off  frequency  of  the  hlter,  0  <  a  <  1.  With  each  frequency  bin  we  associate 
two  filters  with  different  a  values.  As  seen  from  the  second  line  of  Eqn.  (1),  the 
slow  hlter,  with  larger  a  value,  averages  over  more  previous  inputs  x(ri)  than  the  fast 
hlter  with  the  smaller  a  value.  The  slow  hlter  output  is  referred  to  as  the  long-term 
average,  and  the  fast  hlter  output  as  the  short-term  average.  The  slow  hlter  output 
will  track  input  signals  such  as  background  noise  that  changes  slowly  over  time,  and 
reject  shorter  period  signals.  The  fast  hlter  output  will  respond  more  quickly  to  a 
shorter  period  target  signal.  Comparing  the  two  hlter  outputs  in  a  ratio  of  the  fast  to 
slow  will  give  a  measure  that  can  be  used  to  detect  a  target  passage,  while  adjusting 
to  gradual  changes  in  the  background  ambient  noise. 

Both  the  time-domain  and  frequency-domain  Acoustic  Sentinel  algorithms  are  avail¬ 
able  for  acoustic  target  detection  on  the  Starhsh  sensor  nodes.  This  allows  us  to  exper¬ 
imentally  compare  the  detection  performance  of  the  algorithms  to  different  acoustic 
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targets  and  in  different  environmental  conditions.  The  Starfish  sensor  nodes  were 
designed  to  detect  acoustic  targets  in  water  depths  up  to  500  m.  The  frequency  band 
for  detection  in  the  Starfish  sensor  node  is  50  Hz  to  500  Hz.  In  this  band  the  typical 
surface  and  sub-surface  targets  that  the  Starfish  sensor  nodes  were  designed  to  de¬ 
tect  have  characteristically  broad-band  acoustic  signatures.  The  frequency-domain 
Acoustic  Sentinel  was  designed  to  take  advantage  of  the  broad-band  nature  of  the 
expected  target  signatures,  and  reject  any  narrow-band  acoustic  energy  in  the  en¬ 
vironment  that  may  be  present.  A  possible  disadvantage  of  the  frequency-domain 
algorithm  is  that  a  strong  target  with  a  narrow-band  acoustic  spectrum  could  be 
missed,  whereas  the  time-domain  algorithm  essentially  integrates  the  energy  in  the 
whole  band.  Since  both  the  frequency-domain  and  time-domain  Sentinel  algorithms 
can  be  running  simultaneously,  the  detection  differences  between  the  algorithms  could 
be  used  for  possible  target  type  discrimination,  or  to  flag  an  unusual  or  unexpected 
acoustic  environment.  The  frequency-domain  algorithm  allows  for  better  optimiza¬ 
tion  of  detection  performance  because  there  are  more  algorithm  parameters  that  can 
be  tuned  for  particular  target  and  signature  type. 

In  the  frequency-domain  Acoustic  Sentinel  algorithm  we  often  use  1024-point  FFTs 
with  50%  overlap  to  compute  the  average  spectrum.  Typically,  a  total  of  6  FFTs  are 
used  for  the  averaging,  with  the  last  2  FFTs  being  used  in  the  next  group  of  6  FFTs. 
This  results  in  an  average  spectral  power  value  input  to  the  parallel  filters,  in  each 
frequency  bin,  once  every  2048  samples,  or  at  an  effective  sample  rate  of  fs  =  1.22  Hz. 
For  the  slow  and  fast  response  filters,  the  parameter  a  has  the  default  values  of  0.999 
and  0.9,  respectively.  The  frequency  response  of  the  filters  with  these  two  a  values 
is  shown  in  Fig.  2.  The  rise-time  of  the  slow  filter  is  ~  30  minutes,  and  for  the  fast 
filter  ~  20  seconds,  when  data  is  sent  to  the  filters  at  a  rate  of  fs  =  1.22  Hz. 

The  frequency-domain  Acoustic  Sentinel  algorithm  will  trigger  acoustic  target  detec¬ 
tion  when  a  programmable  specified  percentage  (e.g.,  60%)  of  the  frequency  bins  in 
the  50 Hz  to  500 Hz  band  exceed  a  programmable  ratio  threshold  value  (e.g.,  3).  We 
expect  the  frequency-domain  algorithm  to  detect  a  target  if  the  broad-band  signal  to 
noise  ratio  (SNR)  ss  6  dB  or  greater.  A  target  lost  message  is  triggered  when  another 
percentage  (e.g.,  80%)  of  the  bins  fall  below  another  ratio  threshold  parameter  value 
(e.g.,  2).  Using  percentages  of  the  bins  as  target  detect  and  lost  threshold  parameters 
means  that  not  all  of  the  same  bins  have  to  hold  the  required  ratio  threshold  level. 
Some  bins  at  a  certain  time  step  may  fall  below  the  ratio  threshold  value  when  going 
to  the  next  time  step,  but  other  bins  may  increase  over  the  ratio  threshold,  so  that 
the  percentage  of  bins  above  the  target  detect  threshold  is  the  same. 

The  FFT  lengths  and  overlap  values,  the  percentage  of  the  bins  above  and  below 
the  threshold  values,  the  threshold  levels,  and  the  frequency  band  range,  are  all 
algorithm  parameters  that  are  programmable  in  the  Starfish  sensor  node.  Fig.  3 
shows  an  example  of  the  filter  responses  and  resulting  ratio  of  the  fast  to  slow  filter 
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Figure  2:  Normalized  frequency  response  of  the  two  exponential  digital  filters  used 
in  the  Starfish  Acoustic  Sentinel  algorithm. 


outputs  for  a  real  target  pass,  for  one  FFT  bin  at  300  Hz.  This  was  a  particularly 
slow  target  pass  of  a  distant  target,  and  the  filter  slow  and  fast  a  parameters  were 
adjusted  to  0.999  and  0.99  respectively  for  this  target. 

3  Phase  Gradient  bearing  estimation 
algorithm 


The  Phase  Gradient  algorithm  was  originally  designed  as  part  of  a  homing  system 
for  an  Autonomous  Underwater  Vehicle  (AUV)  [2],  In  the  flooded  nose  cone  section 
of  the  AUV  a  tri-axis  cross-dipole,  seven-element  hydrophone  array  was  mounted 
and  used  to  estimate  the  bearing  and  elevation  angles  to  a  single  frequency  beacon 
signal.  A  frequency-domain  bearing  estimation  algorithm  was  the  natural  choice 
to  determine  the  Direction  of  Arrival  (DOA),  to  take  advantage  of  the  increase  in 
SNR  in  going  to  the  frequency-domain  when  processing  a  tonal  signal.  The  on-board 
processing  capability  of  the  AUV  was  unfortunately  not  sufficient  to  implement  the 
Conventional  (Bartlett)  and  Minimum  Variance  Distortionless  Response  (MVDR) 
DOA  algorithms.  This  motivated  the  development  of  the  Phase  Gradient  Algorithm 
for  the  single  frequency  case. 

Both  the  Conventional  and  MVDR  beamforming  algorithms  require  forming  beams 
and  calculating  a  beampower  for  a  large  set  of  discreet  look  directions,  typically 
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Figure  3:  Responses  and  resulting  fast  to  slow  ratio  of  the  two  exponential  digital 
filter  outputs  to  a  target  pass,  for  one  FFT  bin.  The  input  to  the  filters  is  the  average 
spectral  power  for  that  FFT  bin  at  each  time  window. 


searching  over  360°  in  horizontal  bearing  and  from  —90°  to  90°  in  elevation  angle. 
After  calculating  beampowers  over  all  combinations  of  bearing  and  elevation  angles, 
the  single  look  direction  that  gives  the  maximum  output  beampower  with  its  associ¬ 
ated  bearing  and  elevation  angle  define  the  DOA  of  the  incoming  plane  wave.  This 
searching  over  bearing  and  elevation  used  by  the  Conventional  and  MVDR  algorithms 
requires  significant  processing  power.  In  contrast,  the  Phase  Gradient  algorithm  di¬ 
rectly  estimates  the  three  Cartesian  components  of  the  incoming  signal  wavevector 
from  estimated  cross-spectra  between  the  hydrophone  elements.  It  is  robust  against 
hydrophone  failure,  and  every  hydrophone  in  the  array  is  used  to  estimate  each  com¬ 
ponent  of  the  wavevector.  The  Phase  Gradient  algorithm  does  not  require  searching 
over  bearing  and  elevation  angles  and  calculating  a  beampower  at  each  look  direc¬ 
tion,  so  it  very  computationally  efficient.  The  Phase  Gradient  algorithm  was  able  to 
calculate  the  bearing  and  elevation  angles  to  the  homing  beacon  signal  at  the  update 
rate  required  by  the  AUV  homing  system.  This  algorithm  is  well  suited  for  bearing 
estimation  in  embedded  systems  that  have  limited  processing  power. 
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3.1  Phase  Gradient  development 


We  consider  a  planar  acoustic  wave  incoming  to  a  hydrophone  array  from  a  fixed 
direction  of  arrival,  as  shown  in  Fig.  4.  We  define  s(t)  as  the  time  domain  signal 
present  at  hydrophone  #0  due  to  the  plane  wave  arrival,  in  the  absence  of  any  noise. 
We  define  the  vector  ks  as  the  unit  vector  that  points  in  the  direction  of  travel  of  the 
incoming  wave.  We  define  aq(f),  i  =  0 . . .  6  as  the  time  domain  signals  measured  by 
the  hydrophones,  including  noise.  The  time  domain  signal  Xi(t)  measured  at  the  ith 
hydrophone,  in  terms  of  s(t)  is  given  by 


xz(t) 


(2) 


Here  we  have  defined  f%  as  the  position  vector  to  the  ith  hydrophone,  and  c  as  the 
speed  of  sound  in  water.  The  term  rii(t)  represents  the  noise  present  on  the  ith 
hydrophone  that  is  uncorrelated  to  the  incoming  signal.  The  measured  hydrophone 
time  series  Xi(t)  are  delayed  or  advanced  versions  of  s(t)  with  noise  added,  the  delay 
times  depending  on  the  direction  of  arrival  of  the  incoming  signal. 


Figure  4:  The  hydrophone  array  (of  diameter  d )  for  which  the  Phase  Gradient  algo¬ 
rithm  was  designed. 


If  we  take  the  Fourier  Transform  of  Eqn.  (2),  we  obtain  the  frequency  domain  repre¬ 
sentation  of  the  hydrophone  outputs,  given  by 

Xi(f)  =  S{f )  exp  (~jks  o  r;)  +  Nz(f) 

=  \S(f)\exp(-jkaori  +  jlS(f))+Ni(f).  (3) 


6 


DRDC-RDDC-201 5-R025 


The  frequency  domain  representation  of  the  time  signals  is  denoted  by  capital  letters, 
and  j  is  the  imaginary  unit.  The  Argument  of  S(f)  is  denoted  by  IS(f).  The 
incoming  signal  wavevector  ks  =  \ks\ks  =  — —ks  of  Eqn.  (3)  is  a  function  of  frequency. 
From  Eqn.  (3)  we  can  see  that  for  a  given  frequency,  the  spatial  phase  at  any  position 
vector  r  of  the  incoming  signal,  neglecting  noise,  is  given  by 

ip(r)  =  -ka  o  f  +  IS(f)  .  (4) 

The  phase  angle  IS(f)  does  not  depend  upon  the  spatial  position  vector  r,  and  will 
have  zero  gradient.  If  we  can  estimate  the  gradient  of  the  spatial  phase  "0(A)  we  will 
have  an  estimate  of  the  wavevector  V0(r)  =  —  ks,  and  thus  have  the  arrival  direction. 


3.2  Orthogonal  array  Phase  Gradient 


We  estimate  the  gradient  of  the  spatial  phase  given  by  Eqn.  (4),  for  the  cross-dipole 
array  of  Fig.  4,  from  the  Argument  of  a  modified  definition  of  the  cross-spectrum 
between  two  hydrophones.  We  define  Xt  as  being  the  complex  FFT  spectral  value 
in  a  particular  frequency  bin  for  the  ith  hydrophone.  Here  we  write  ks  in  Cartesian 
coordinates  ks  =  (kx,ky,kz),  use  (  )  to  denote  expectation  value,  and  0,  to  denote 
the  phase  of  the  ith  hydrophone.  With  the  array  diameter  being  denoted  by  d,  we 
have  for  the  rc-axis  direction  a  "new"  cross-spectrum  defined  as 


(xiX*2)new  =  (X1{X0  +  X3  +  X4  +  X5  +  X6)*  +  (X0  +  X3  +  X4  +  X3  +  X6)X*) 
=  2(SS*)  ^1  +  2  cos^ky^)  +  2  cos(/c2^)^  eJ'5 Wa-VA  _  (5) 


In  Eqn.  (5)  the  phase  contributions  from  the  y  and  £-axis  directions  combine  to 
create  a  real-valued  amplitude  factor  depending  on  ky  and  kz  because  of  the  array 
symmetry,  leaving  only  half  the  phase  difference  between  hydrophones  #1  and  ^2  as 
the  Argument  of  the  complex  exponential.  Using  the  "new"  cross-spectrum  ensures 
that  all  the  array  hydrophones  are  used  in  the  estimate  of  the  phase  difference.  We 
can  estimate  the  ^-component  of  V0(r)  from 

=  v-i-fe  =  a±  ,, .. 

d  d  dx 

For  the  y  and  £-axis  directions  we  can  perform  similar  operations  to  Eqs.  (5)  and  (6), 
obtaining  the  remaining  components  of  V0(r)  =  —  ks  and  the  direction  of  arrival.  We 
estimate  the  bearing  and  elevation  angles  from  the  wavevector  components  in  each 
frequency  bin  from 


9  =  arctan 


(j)  =  arctan 


(7) 
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3.3  Tetra  array  Phase  Gradient 

The  Phase  Gradient  algorithm  has  been  adapted  for  use  with  a  4-element  tetrahedral 
array  shape  [3].  Fig.  5  shows  the  Tetra  array  geometry,  where  d  =  0.162  m  is  the 
distance  from  the  origin  to  each  x  —  y  plane  hydrophone,  and  dz  =  0.228  m  is  the 
height  of  the  top  hydrophone  above  the  x  —  y  plane.  We  estimate  the  arrival  angles 
from  the  components  of  the  wavevector  using  Eq.  (7).  The  wavevector  components 
for  the  Tetra  geometry  are  given  by  Eq.  (8). 


Figure  5:  The  Tetra  hydrophone  array  for  which  the  Phase  Gradient  algorithm  was 
adapted. 


kx  =  ^Arg{(XoX^(X0X*)^{X3X^{X3X:)^{X3X*)^}, 


k,,  = 


-^Arg  {(XoXD^XzX^iXzX^iXs X*)Hx2X*3)*} , 


(8) 


kz  =  -Arg{(X0X*^(X0X*2^(X0X*3)^}. 
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4  Comparison  of  Phase  Gradient  to  other 
DOA  algorithms 


When  considering  DOA  algorithms  to  implement  on  the  second  generation  Starfish 
sensor  nodes  the  Phase  Gradient  algorithm  and  MVDR  algorithms  were  chosen.  The 
Phase  Gradient  algorithm  was  chosen  because  of  its  good  performance  in  the  AUV 
homing  system  at  a  single  frequency.  This  was  a  good  opportunity  to  extend  the  use  of 
the  algorithm  to  broad-band  signals,  and  to  experimentally  evaluate  its  performance 
on  a  7-element  cross-dipole  array  shape.  The  MVDR  algorithm  was  chosen  because 
of  its  performance  in  simulations,  and  its  ability  to  resolve  more  than  one  DOA  signal 
simultaneously.  Both  algorithms  are  similar  in  that  they  both  are  frequency-domain 
algorithms  that  result  in  a  DOA  bearing  estimate  in  each  FFT  frequency  bin.  In 
computer  simulations  the  performance  of  both  of  these  algorithms  was  compared  to 
the  performance  of  the  standard  time-domain  crossed-dipole  DOA  algorithm  used  in 
DIFAR  sonobuoys.  The  DIFAR  time-domain  algorithm  [4]  will  be  referred  to  here  as 
the  DIFAR  TD  algorithm. 

In  the  computer  simulations  the  dipole  array  diameter  was  set  to  0.4  m,  the  sound 
speed  1490  m/s,  the  sample  rate  8000  samples  per  second.  Time  blocks  of  10  second 
duration  were  used,  and  for  the  frequency-domain  algorithms  the  FFT  lengths  were 
set  to  8000  points.  The  required  average  spectra  were  calculated  using  50%  overlap 
for  a  total  of  19  averages  per  time  block.  The  plane  wave  arrival  signal  and  the  noise 
on  each  hydrophone  were  simulated  as  broad-band  gaussian  white  noise  sequences. 
A  single  10  second  time  block  with  80000  samples  is  considered  as  one  simulation, 
which  produces  one  bearing  estimate  (only  horizontal  bearings  are  considered  in  this 
paper)  from  each  DOA  algorithm.  The  true  DOA  bearing  for  the  simulations  was  set 
to  45°.  All  the  processing  was  performed  in  the  frequency  band  of  50  Hz  to  500  Hz. 

4.1  Average  Histogram  method 

For  both  the  Phase  Gradient  and  MVDR  algorithms  we  get  an  estimate  of  the  DOA 
of  the  incoming  plane  wave  signal  at  each  FFT  frequency  bin  in  the  band  of  interest. 
Two  approaches  were  used  to  estimate  a  single  DOA  value  from  the  frequency  domain 
information  in  all  the  bins. 

The  first  approach  was  to  calculate  the  bearing  in  each  frequency  bin  and  then  create  a 
histogram  of  the  bearings.  We  divide  the  horizontal  bearing  range  of  —180°  to  +180° 
into  a  set  of  equally  spaced  bearing  bins  of  1°  width.  For  one  simulation  we  get  a 
bearing  estimate  at  each  FFT  frequency  bin.  For  every  frequency  bin,  we  add  one 
count  to  the  bearing  bin  whose  associated  bearing  range  contains  the  DOA  estimate. 
The  sum  of  all  the  counts  over  all  the  bearing  bins  equals  the  number  of  FFT  frequency 
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Degrees 

Figure  6:  Phase  Gradient  bearing  histogram,  un-smoothed  (top)  and  smoothed  (bot¬ 
tom)  for  one  simulation  at  12  dB  SNR  and  DO  A  of  45°.  A  Gaussian  shape  is  fitted 
to  the  histograms. 


bins.  Fig.  6  shows  a  Phase  Gradient  bearing  histogram,  un-smoothed  and  smoothed 
for  a  single  simulation,  with  the  true  DOA  set  at  45°  and  the  SNR  12  dB.  We  can 
see  that  the  smoothed  histogram  is  approaching  a  Gaussian  distribution  shape.  The 
arctangent  DOA  estimator  is  an  unbiased  estimator  [5],  and  its  probability  density 
function  (pdf)  is  estimated  from  the  smoothed  histogram.  A  21-point  boxcar  average 
was  used  for  smoothing,  which  creates  real-valued  histogram  counts.  To  estimate 
the  DOA  we  then  perform  a  sum  of  the  pdf  times  the  bearing  value  over  a  31°  width 
about  the  peak.  Summing  over  the  peak  prevents  bearings  that  are  too  far  away  from 
the  peak  to  be  included  in  the  DOA  estimate.  The  same  procedure  is  performed  for 
the  MVDR  algorithm  with  similar  results,  although  the  Phase  Gradient  algorithm 
estimated  pdf  has  a  noticeably  smaller  variance.  This  method  of  DOA  estimation 
will  be  referred  to  as  the  Average  Histogram  (Ave  Hist)  method. 
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4.2  Average  Direction  Cosine  method 

The  second  method  for  estimating  the  DOA  for  the  Phase  Gradient  and  MVDR 
algorithms  was  to  calculate  the  average  direction  cosine  values  to  each  Cartesian  axis 
over  all  the  frequency  bins,  and  then  perform  the  required  arctangent  operation.  For 
horizontal  bearings  considered  here,  this  is  just  averaging  over  the  frequency  bins  the 
sin(0)  and  cos(0)  values  of  the  DOA  angle  9  estimated  in  each  frequency  bin.  The 
arctangent  is  then  computed  using  these  averaged  values.  This  procedure  will  be 
referred  to  as  the  Average  Direction  Cosine  (Ave  DC)  method. 

4.3  Simulation  results 

To  estimate  the  statistics  of  the  different  DOA  algorithms  we  performed  100  simu¬ 
lations  on  each  algorithm,  each  with  different  seed  values  for  the  Gaussian  sequence 
generators,  for  every  SNR  value  in  the  range  from  -9  dB  to  +12  dB  in  3  dB  steps. 
Fig.  7  shows  the  resulting  DOA  estimates  from  all  the  algorithms  for  100  simulations, 
with  SNR  -3  dB  and  the  true  DOA  at  45°.  Fig.  8  shows  the  estimated  DOA  mean 
values  and  corresponding  error  bars  of  ±1  standard  deviation,  for  all  algorithms,  as 
a  function  of  SNR. 
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Figure  7:  Estimates  of  the  DOA  from  all  algorithms  for  100  simulations,  with  SNR 
—3  dB  and  the  true  DOA  at  45°. 

We  can  see  from  Fig.  8  that  apart  from  the  MVDR  Ave  Hist  algorithm,  the  rest  of  the 
DOA  algorithms  give  similar  results  for  SNR  values  above  0  dB.  Above  12  dB  SNR 
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Figure  8:  Mean  and  Standard  Deviations  of  the  estimated  DO  A  for  100  simulations, 
as  a  function  of  Signal  to  Noise  Ratio  (SNR).  The  true  DO  A  is  at  45°,  indicated  by 
the  dashed  line. 

we  have  the  standard  deviation  of  all  the  DOA  algorithm  estimates  at  or  below  ~  1°. 
As  the  SNR  value  decreases  below  0  dB  we  can  see  that  the  standard  deviations  of 
the  DOA  estimates  increase  rapidly.  The  Phase  Gradient  Ave  DC  algorithm  and  to  a 
slightly  lesser  extent  the  Phase  Gradient  Ave  Hist  algorithm  have  similar  performance 
to  the  time-domain  DIFAR  TD  method,  showing  similar  means  and  standard  devia¬ 
tions  over  the  whole  range  of  SNR  values.  For  SNR  values  below  0  dB  the  MVDR  Ave 
DC  algorithm  consistently  has  the  smaller  standard  deviation  of  all  the  algorithms. 

5  Phase  Gradient  Histogram  Tracking 
example 


Figs.  9  and  10  show  the  results  of  the  Phase  Gradient  processing  and  Histogram 
Tracking  using  data  collected  from  a  Tetra  array.  The  data  was  supplied  by  the 
German  Federal  Defence  Organization  Bundeswehr,  under  the  Next  Generation  Au¬ 
tonomous  Systems  (NGAS)  NATO  Joint  Research  Project  (JRP)  agreement.  The 
data  was  collected  at  the  NGAS  2012  engineering  trial  in  Horten,  Norway.  The  data 
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Figure  9:  Horizontal  bearings  from  the  Phase  Gradient  algorithm,  with  color  coded 
bearing  as  a  function  of  time  and  frequency. 


set  contains  a  broad-band  acoustic  target  pass.  For  this  data  set  the  hydrophone  sam¬ 
ple  rate  is  20000  Hz.  Data  was  processed  from  50  Hz  to  2600  Hz  using  10000-point 
FFTs.  The  frequency  band  upper  limit  of  the  processing  is  determined  by  the  Tetra 
array  size.  The  Phase  Gradient  algorithm  requires  the  maximum  array  element  sep¬ 
aration  to  be  less  than  half  a  wavelength.  The  horizontal  bearing  range  was  divided 
into  360  bearing  bins,  each  1°  wide.  A  total  of  15  FFTs  were  used  for  the  spectral 
averaging  with  no  overlap,  resulting  in  480  time  windows  each  7.5  seconds  long.  In 
Fig.  9  we  see  the  Phase  Gradient  results  as  color  coded  horizontal  bearing  as  a  func¬ 
tion  of  time  and  frequency.  Fig.  10  shows  the  Histogram  Tracking  results,  with  color 
coded  histogram  counts.  The  broad-band  target  track  for  this  data  set  is  well  defined 
using  the  Phase  Gradient  algorithm  and  Histogram  Tracking  procedure,  and  either 
the  Average  Direction  Cosine  or  Average  Histogram  method  gives  a  good  single  value 
estimate  for  the  DOA.  Determining  when  to  run  the  DOA  estimate  algorithms  could 
be  signalled  from  the  Sentinel  algorithms  once  a  target  is  detected,  or  a  threshold 
level  on  the  histogram  count  could  also  be  used. 
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Figure  10:  Histogram  tracking  of  the  horizontal  bearing  using  the  Phase  Gradient 
algorithm,  with  color  coded  histogram  counts. 

6  Conclusion 


On  the  DRDC  Starfish  sensor  nodes  we  have  time-domain  and  frequency-domain  Sen¬ 
tinel  algorithms  for  target  detection  implemented.  The  frequency-domain  Acoustic 
Sentinel  algorithm  was  designed  to  detect  broad-band  acoustic  targets.  The  tirnedo- 
main  Acoustic  Sentinel  algorithm  integrates  all  energy  in  the  band  of  interest  and 
does  not  depend  on  the  spectral  shape  of  the  signal.  Since  both  can  be  running 
simultaneously,  differences  in  their  detections  could  be  used  to  flag  different  target 
types  or  unusual  acoustic  environments. 

The  Phase  Gradient  algorithm  performance  in  determining  the  DOA  of  a  broad-band 
signal  was  evaluated  using  computer  simulations  and  tested  experimentally.  Two 
methods  were  examined  that  combined  the  bearing  information  from  each  FFT  fre¬ 
quency  bin  into  one  DOA  estimate;  the  Average  Direction  Cosine  (Ave  DC)  and  the 
Average  Histogram  (Ave  Hist)  methods.  Results  from  the  Phase  Gradient  algorithm 
were  compared  to  those  from  the  MVDR  and  the  standard  time-domain  DIFAR  TD 
DOA  algorithms.  The  Phase  Gradient  algorithm  using  the  Ave  DC  method  gave  simi¬ 
lar  DOA  estimates  with  similar  standard  deviations  as  the  DIFAR  TD  algorithm  over 
all  SNR  values.  The  Phase  Gradient  Ave  Hist  algorithm  gave  slightly  larger  standard 
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deviations  below  0  dB  than  the  DIFAR  TD.  Above  0  dB  SNR  all  the  algorithms 
performed  well.  Below  0  dB  SNR  the  standard  deviations  of  the  DOA  estimates 
increased,  with  the  MVDR  Ave  DC  algorithm  having  the  smallest  standard  devia¬ 
tion,  while  the  other  algorithms  performed  similarly,  except  for  the  MVDR  Ave  Hist 
algorithm  which  had  a  consistently  larger  standard  deviation  for  reasons  unknown. 

For  detecting  a  broad-band  DOA  signal  the  DIFAR  TD,  Phase  Gradient  Ave  DC, 
and  MVDR  Ave  DC  algorithms  all  performed  similarly  above  0  dB,  with  the  MVDR 
Ave  DC  algorithm  having  smaller  standard  deviations  below  0  dB.  For  an  embedded 
system  with  limited  processing  power  the  Phase  Gradient  Ave  DC  and  DIFAR  TD 
algorithms  would  have  the  advantage  over  the  MVDR  algorithm  because  they  require 
much  less  processing  power.  Since  the  Phase  Gradient  algorithm  is  applied  in  the 
frequency  domain,  it  would  have  the  advantage  if  processing  of  particular  frequency 
bands  was  required.  One  such  example  might  be  the  tracking  of  a  target  using  its 
tonal  lines  and  associated  harmonics. 

The  Phase  Gradient  algorithm  has  the  disadvantage  of  only  being  applicable  to  very 
symmetric  array  shapes,  such  as  the  crossed-dipole  or  tetrahedral  arrays.  The  MVDR 
algorithm  can  be  applied  to  a  general  array  shape,  and  also  has  the  advantage  of 
being  able  to  resolve  the  DOA  of  more  than  one  arrival  simultaneously.  The  Phase 
Gradient  algorithm  is  much  more  efficient  than  the  MVDR  algorithm  and  is  more 
applicable  to  embedded  systems  with  limited  processing  power.  In  a  tracking  system 
the  Phase  Gradient  algorithm  would  be  simple  to  implement  because  it  automatically 
gives  the  bearing  without  searches.  The  MVDR  algorithm  may  require  a  complex 
tracking  algorithm  to  follow  the  target  from  update  to  update  to  reduce  processing 
requirements  by  only  searching  over  look  directions  near  the  DOA  from  the  previous 
update. 

The  frequency-domain  and  histogram  tracking  results  presented  in  Figs.  9  and  10  from 
the  Phase  Gradient  algorithm  give  an  experienced  observer  a  great  deal  of  information 
about  the  frequency  content  and  tracking  effectiveness  for  a  target.  These  would 
be  useful  real-time  displays  when  tracking  a  target  from  a  shore  station  where  the 
underwater  sensors  would  be  cabled  to  shore.  Because  of  the  efficiency  of  the  Phase 
Gradient  algorithm  it  may  be  more  possible  to  create  real-time  updated  displays  such 
as  these  on  typical  computers  than  with  the  MVDR  algorithm. 


DRDC-RDDC-201 5-R025 


15 


This  page  intentionally  left  blank. 


16 


DRDC-RDDC-201 5-R025 


References 


[1]  Otnes,  R.  (May  2012),  NILUS  -  An  Underwater  Acoustic  Network 
Demonstration  System,  10th  International  Mine  Warfare  Technology 
Symposium.  Monterey,  California,  USA. 

[2]  Lucas,  C.,  Heard,  G.,  and  Pelavas,  N.  (2013),  The  Phase  Gradient  Bearing 
Estimation  Algorithm,  Proceedings  of  Meetings  on  Acoustics  (POMA). 

[3]  Datta,  U.,  Otnes,  R.,  and  Lucas,  C.  (Sep.  2010),  Bearing  estimation  using  small 
tetrahedral  passive  hydrophone  array,  Oceans  Conference.  Seattle  Wash. 

[4]  Davies,  S.  (Sep.  1987),  Bearing  accuracies  for  arctan  processing  of  crossed  dipole 
arrays,  Oceans  Conference.  Halifax,  NS. 

[5]  Maranda,  B.  (Sep.  2003),  The  statistical  accuracy  of  an  arctangent  estimator, 
Oceans  Conference.  Pages(s)  2127-2132  Vol.4. 


DRDC-RDDC-201 5-R025 


17 


This  page  intentionally  left  blank. 


18 


DRDC-RDDC-201 5-R025 


DOCUMENT  CONTROL  DATA 

(Security  markings  for  the  title,  abstract  and  indexing  annotation  must  be  entered  when  the  document  is  Classified  or  Designated.) 


1 .  ORIGINATOR  (The  name  and  address  of  the  organization  preparing  the 
document.  Organizations  for  whom  the  document  was  prepared,  e.g.  Centre 
sponsoring  a  contractor’s  report,  or  tasking  agency,  are  entered  in  section  8.) 

2a.  SECURITY  MARKING  (Overall  security 
marking  of  the  document,  including 
supplemental  markings  if  applicable.) 

DRDC  -  Atlantic  Research  Centre 

PO  Box  1012,  Dartmouth  NS  B2Y  3Z7,  Canada 

UNCLASSIFIED 

2b.  CONTROLLED  GOODS 

(NON-CONTROLLED 

GOODS) 

DMC  A 

REVIEW:  GCEC  APRIL  2011 

3.  TITLE  (The  complete  document  title  as  indicated  on  the  title  page.  Its  classification  should  be  indicated  by  the  appropriate 
abbreviation  (S,  C  or  U)  in  parentheses  after  the  title.) 


DRDC  Starfish  Acoustic  Sentinel  and  Phase  Gradient  Histogram  Tracking 


4.  AUTHORS  (Last  name,  followed  by  initials  -  ranks,  titles,  etc.  not  to  be  used.) 

Lucas,  C.  E.;  Heard,  G.  J.;  Pelavas,  N. 


5.  DATE  OF  PUBLICATION  (Month  and  year  of  publication  of 
document.) 

6a.  NO.  OF  PAGES  (Total 
containing  information. 

Include  Annexes, 

Appendices,  etc.) 

6b.  NO.  OF  REFS  (Total 
cited  in  document.) 

April  2015 

26 

5 

7.  DESCRIPTIVE  NOTES  (The  category  of  the  document,  e.g.  technical  report,  technical  note  or  memorandum.  If  appropriate,  enter 
the  type  of  report,  e.g.  interim,  progress,  summary,  annual  or  final.  Give  the  inclusive  dates  when  a  specific  reporting  period  is 
covered.) 


Scientific  Report 


8.  SPONSORING  ACTIVITY  (The  name  of  the  department  project  office  or  laboratory  sponsoring  the  research  and  development  - 
include  address.) 

DRDC  -  Atlantic  Research  Centre 

PO  Box  1012,  Dartmouth  NS  B2Y  3Z7,  Canada 


9a.  PROJECT  OR  GRANT  NO.  (If  appropriate,  the  applicable 
research  and  development  project  or  grant  number  under 
which  the  document  was  written.  Please  specify  whether 
project  or  grant.) 

9b.  CONTRACT  NO.  (If  appropriate,  the  applicable  number  under 
which  the  document  was  written.) 

01CL01 

10a.  ORIGINATOR’S  DOCUMENT  NUMBER  (The  official 

document  number  by  which  the  document  is  identified  by  the 
originating  activity.  This  number  must  be  unique  to  this 
document.) 

tOb.  OTHER  DOCUMENT  NO(s).  (Any  other  numbers  which  may 
be  assigned  this  document  either  by  the  originator  or  by  the 
sponsor.) 

DRDC-RDDC-201 5-R025 

1 1 .  DOCUMENT  AVAILABILITY  (Any  limitations  on  further  dissemination  of  the  document,  other  than  those  imposed  by  security 
classification.) 


(X)  Unlimited  distribution 

(  )  Defence  departments  and  defence  contractors;  further  distribution  only  as  approved 
(  )  Defence  departments  and  Canadian  defence  contractors;  further  distribution  only  as  approved 
(  )  Government  departments  and  agencies;  further  distribution  only  as  approved 
(  )  Defence  departments;  further  distribution  only  as  approved 
(  )  Other  (please  specify): 


12.  DOCUMENT  ANNOUNCEMENT  (Any  limitation  to  the  bibliographic  announcement  of  this  document.  This  will  normally  correspond 
to  the  Document  Availability  (11).  However,  where  further  distribution  (beyond  the  audience  specified  in  (11))  is  possible,  a  wider 
announcement  audience  may  be  selected.) 


1 3.  ABSTRACT  (A  brief  and  factual  summary  of  the  document.  It  may  also  appear  elsewhere  in  the  body  of  the  document  itself.  It  is  highly 
desirable  that  the  abstract  of  classified  documents  be  unclassified.  Each  paragraph  of  the  abstract  shall  begin  with  an  indication  of  the 
security  classification  of  the  information  in  the  paragraph  (unless  the  document  itself  is  unclassified)  represented  as  (S),  (C),  or  (U).  It  is 
not  necessary  to  include  here  abstracts  in  both  official  languages  unless  the  text  is  bilingual.) 

DRDC  -  Atlantic  Research  Centre,  an  agency  of  the  Canadian  Department  of  National  Defence, 
has  developed  a  prototype  multi-influence  underwater  surveillance  sensor  called  the  Starfish. 
The  Starfish  has  acoustic,  electric,  and  magnetic  field  sensors,  and  onboard  signal  processing 
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the  array  centre. 

Acoustic  target  detection  is  performed  using  both  time-domain  and  frequency-domain  Sentinel 
algorithms.  The  Sentinel  algorithms  use  parallel  exponential  filters,  with  the  frequency-domain 
algorithm  using  parallel  filters  in  each  frequency  bin. 

A  Phase  Gradient  bearing  estimation  algorithm  is  used  to  perform  the  acoustic  target  tracking. 
This  is  a  frequency-domain  algorithm  that  estimates  the  gradient  of  the  spatial  phase  of  the 
incoming  acoustic  signal,  and  directly  computes  the  wavevector  of  the  acoustic  arrival  in  each 
frequency  bin.  A  Histogram  Tracking  procedure  is  used  to  display  and  determine  the  track  angles 
for  a  broad-band  acoustic  target.  This  paper  will  discuss  the  details  of  the  frequency-domain 
Acoustic  Sentinel  algorithm  and  the  Phase  Gradient  bearing  estimation  algorithm  with  Histogram 
Tracking. 
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